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Abstract 

Data-dependent  interpolatory  techniques  can  be  used  in  the  reconstruction  step  of  a 
multiresolution  “a  la  Marten" . These  interpolatory  techniques  lead  to  nonlinear  mul- 
tiresolution schemes.  When  dealing  with  nonlinear  algorithms,  the  issue  of  the  stability 
needs  to  be  carefully  considered.  In  this  paper  we  analyze  and  compare  several  strategies 
for  image  compression  and  their  ability  to  effectively  control  the  global  error  due  to 
compression. 


1 Introduction 

Multiscale  transformations  are  being  used  in  recent  times  in  the  first  step  of  transform 
coding  algorithms  for  image  compression.  Ideally,  a multiscale  transformation  allows  for 
an  efficient  representation  of  the  image  data,  which  is  then  processed  using  a (non- 
reversible)  quantizer  and  passed  on  to  the  encoder  which  produces  the  final  compressed 
set  of  data  which  is  ready  to  be  transmitted  or  stored.  Compression  is  indeed  achieved 
during  the  second  and  third  steps:  the  quantization  and  the  encoding  of  the  transformed 
set  of  discrete  data. 

It  is  quite  clear  that  the  properties  of  the  multiscale  transformation  are  most  im- 
portant in  the  overall  performance  of  the  transform  coding  algorithm.  Until  recently,  the 
multiscale  transformations  used  for  image  compression  were  always  based  on  linear  filter 
banks,  however,  the  nonlinear  alternative  has  been  explored  lately  by  various  authors 
from  different  points  of  view,  and  preliminary  results  show  the  alternative  to  be  very 
promising  [12,  8,  6,  2,  3].  The  key  question  when  using,  or  even  designing,  a nonlinear 
multiscale  transformation  is  that  of  stability.  In  order  for  such  transformations  to  be 
useful  tools  in  image  coding,  it  is  absolutely  necessary  to  keep  a tight  control  on  the 
effect  of  quantization  errors  in  the  decoding  process. 

In  this  paper  we  examine  the  question  of  stability  for  nonlinear  multiscale  trans- 
formations within  Harten’s  framework  for  multiresolution  [14,  15].  Harten’s  framework 
is  broad  enough  to  include  all  classical  wavelet  transformations  as  particular  cases  (just 
as  it  happens  in  the  Lifting  framework  of  W.  Sweldens  [17],  developed  slightly  later  in 
time  but  independently),  however  the  design  of  the  multiscale  transformation  is  done 
directly  on  the  spatial  domain. 
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The  building  blocks  of  Harten’s  multiresolution  framework  are  two  operators  that 
connect  adjacent  resolution  levels.  The  Decimation  (or  also,  Restriction ) operator  is  a 
linear  operator  which  acts  as  a low-pass  filter,  extracting  low-resolution  information  from 
a discrete  data  set.  The  Prediction  operator  (also  Projection)  uses  low-resolution  data  to 
predict  discrete  data  at  a higher  resolution  level.  It  is  precisely  the  design  of  this  operator 
what  distinguishes  Harten’s  framework  from  all  other  multiresolution  frameworks.  The 
prediction  operator  is  based  on  a consistent  Reconstruction  technique,  and  this  opens  up 
a tremendous  number  of  possibilities  in  the  design  of  multiresolution  schemes.  The  use 
of  the  reconstruction  process  as  a design  tool  makes  it,  conceptually,  a simple  matter 
to  introduce  adaptivity  into  the  multiscale  transformation;  we  only  need  to  make  the 
reconstruction  process  data-dependent  [5,  4,  14]. 

This  paper  is  organized  as  follows.  In  Section  2 we  recall  the  so-called  cell-average 
framework,  an  appropriate  setting  for  image  compression,  and  describe  a class  of  nonlin- 
ear prediction  operators  obtained  by  mean-average  interpolation  [10, 14, 15].  In  Section  3 
we  examine  the  question  of  stability  for  nonlinear  multiscale  transformations  and  relate 
it  to  the  synchronization  of  the  data-dependent  choices  made  in  the  encoder  and  the 
decoder.  We  also  include  a set  of  numerical  experiments  that  illustrate  he  performance 
of  several  nonlinear  multiscale  transformations. 


2 Multiscale  transformations  in  the  cell-average  setting 

Harten’s  general  framework  for  multiresolution  [15]  relies  on  two  operators,  Decimation 
and  Prediction,  that  define  the  basic  interscale  relations.  These  operators  act  on  finite 
dimensional  linear  vector  spaces,  Vj,  that  represent  the  different  resolution  levels  (j 
increasing  implies  more  resolution) 

(a)  Dj  : Vj  (b)  Pj  : V^1  — > Vj , (2.1) 


and  must  satisfy  two  requirements  of  algebraic  nature;  D-7  needs  to  be  a linear  operator 
and  D^Pj  = IVj- 1,  i.e.,  the  identity  operator  on  the  lower  resolution  level  represented 
by  Vi-1.  For  all  practical  purposes,  can  be  considered  as  spaces  of  finite  dimensional 
sequences. 

Using  these  two  operators,  a vector  (i.e.,  a discrete  sequence)  u7  G Vi  can  be  decom- 
posed and  reassembled  as  follows 


(a) 


v- 


•i  —>  vi  1 

pj 


Diyi 

vi  - Pjvi-1 


(b)  vi  = Pjyi  1 + e-7  (2.2) 


where  e7  represents  the  error  in  trying  to  predict  the  jth  level  data,  vi,  from  the  low 
resolution  data  v-7-1  = Diyi,  using  the  prediction  operator  Pj  . 

In  the  cell-average  setting,  the  discrete  data  are  interpreted  as  the  cell-averages  of 
a function  on  an  underlying  grid,  which  determines  the  level  of  resolution  of  the  given 
data.  The  one  dimensional  case,  in  which  one  considers  a set  of  nested  dyadic  grids  on 
the  interval  [0, 1],  {X7},  j > 0 of  size  hj  = 2_-7/i0, 


Xi  = {x{}  x{  = i ■ hj,  i = 0, . . . , Nj  Nj  ■ hj  = 1 


(2.3) 
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is  the  easiest  one  to  describe,  and  it  is  also  directly  applicable  to  two-dimensional  (2D) 
data  via  tensor  product  [2,  3]  (the  cell-average  framework  in  several  dimensions  and 
non- tensor  product  (unstructured)  grids  is  considered  in  e.g.  [1]). 

In  this  simple  one-dimensional  setting,  the  cell-average  framework  is  characterized  by 
the  following  decimation  operator 

pV)<  = i(t4_i  + t4i).  (2.4) 

where  Nj  is  the  number  of  equally  spaced  intervals  on  , the  grid  on  [0, 1]  that  repres- 
ents the  jth  resolution  level.  The  consistency  requirement  for  the  prediction  operator, 
i.e.,  D^Pj  = IVj- 1 which  is  the  only  necessary  requirement  for  the  prediction  in  Harten’s 
framework,  becomes  then 

(Py-1)*-!  + {PjVj-l)2i  = 2 vj-\  (2.5) 

Observe  that  (2.4)  and  (2.5)  imply  that 


( V"1)®  + {Pjv’-'hi-i  = 4 + 4-i* 

Hence 

4-1  = 4-i  - (Pjvi-'hi-i  = (/>i_1)2*  - 4 = -4- 


Therefore  the  prediction  errors  at  even  and  odd  grid  points  on  the  jth  level  in  (2.2) 
are  not  independent.  By  considering  only  the  prediction  errors  at  (for  example)  the 
odd  points  of  the  grid  X-\  one  immediately  gets  a one-to-one  correspondence  between 
the  sets  {v*}^  *->  {{vf1}^-1,  with  d\  = e?2i_x  and  = DjvP  The 

one-dimensional  multiscale  transformation  and  its  inverse  can  be  written  as  follows, 


vL  — + MvL  = (v°,d\...,dL) 


vd  = (v°,  dl , . . . , dL)  — + M 1vd 


i 


For  j = L, . . . , 1 
For  i — 1, . . . , Nj- 1 

= (4+4-i)/2 
di  = 4-1  - 

( For  j — 1,. . . ,L 
| For  i = 1, . . . , Nj-i 

1 4-i  = (Pf^'hi-i  + d? 

{ 4,  = 2v]-l-v32i_l 


>-  (2-6) 


)■  (2.7) 


Observe  that  since  dj  = eJ2i_1  = -eJ2l,  the  consistency  relation  (2.5)  implies  that  the 
computation  of  v2i  in  (2.7)  is  equivalent  to 

4 = 2vl~l  - 4-1  = (Pjvi~lhi  -4  = (7>i~1)2t  + 4‘  (2.8) 

Therefore  (2.6)  and  (2.7)  are  just  the  repeated  application  of  the  decomposition  and  reas- 
sembling specified  in  (2.2) (a)  and  (2.2) (b).  Thus  (2.6)  defines  a multiscale  transformation 
and  (2.7)  is  the  inverse  transformation,  whether  or  not  the  prediction  operator  is  linear. 

Next,  we  follow  [4,  14,  15]  to  describe  a class  of  linear  prediction  operators  that  leads 
to  the  (1,  M)  branch  of  the  Cohen-Daubechies-Feauveau  family  [7],  which  is  biorthogonal 
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to  the  box  function  [11,  15].  This  class  is  also  considered  in  [6]  within  the  lifting  frame- 
work, where  it  is  described  as  a particular  case  of  Donoho’s  average  interpolation  [9]. 


Given  an  integer  s > 1,  for  each  1 < * < Nj-i  we  construct  a polynomial,  pt(x),  of 
degree  2s  such  that 


1 

hj~i 


L 


i+l 


1 

'i+l  — 1 


Pi(x)dx 


J-1 

vi+l  > 


for  l = —s, ...  ,s. 


(2.9) 


There  are  various  ways  to  prove  that  Pi(x)  in  (2.9)  always  exists  and  it  is  uniquely 
defined  by  the  2s  + 1 conditions  in  (2.9)  [1,  9,  14].  Then  we  define 


(py-%  = 


(2.10) 


The  prediction  operator  defined  by  (2.10)  is  data-independent,  hence  linear,  and 
it  clearly  satisfies  the  consistency  relation  (2.5).  It  can  be  shown  that  the  multiscale 
transformations  (2.6)  and  (2.7)  for  this  class  of  prediction  operators  turns  out  to  be  the 
(1  ,M  = 2s  + 1)  branch  of  the  Cohen-Daubechies-Feauveau  family. 

A nonlinear  prediction  operator  is  obtained  if  we  construct  Pi(x)  in  a data-dependent 
way.  An  example  of  nonlinear  multiresolution  transformation  constructed  in  this  fashion 
is  considered  in  [14,  4,  2],  where  a nonlinear  ENO-type  technique  (Essentially  Non  Os- 
cillatory, see  [16])  is  used  to  construct  Pi(x).  The  key  idea,  which  is  in  essence  common 
to  the  approach  used  in  designing  nonlinear  filter  banks,  is  to  avoid  using  data  across 
an  edge  for  the  prediction  step. 

The  ENO  nonlinear  technique  is  better  described  if  we  associate  to  each  polynomial 
piece  Pi(x)  a stencil,  Si,  which  is  the  set  of  indices  of  the  values  used  to  define  Pi(x).  In 
the  linear  case  <S,  = {i  — s, . . . ,i  + s};  the  stencil  is  independent  of  the  data  set  {vJ~\} 
and,  as  a consequence,  Pj  is  a linear  operator.  In  the  ENO  technique  described  in  [16],  the 
selection  of  stencil  is  made  in  a data-dependent  way  using  the  divided  differences  of  the 
data  as  a measure  of  its  smoothness.  Large  divided  differences  occur  when  considering 
data  across  an  edge,  while  divided  (or  undivided)  differences  of  data  on  smoother  regions 
tend  to  be  smaller  in  size. 


The  information  contained  in  the  divided  differences  is  then  used  to  decide  what  is 
Si  for  each  i,  with  the  only  restriction  that  i & Si  (to  satisfy  the  consistency  requirement 
(2.5)).  We  follow  [4]  and  consider  all  polynomial  pieces  of  the  same  degree.  In  our  case 
:tf  Sl  = 2s,  but  in  principle  one  could  decide  to  lower  the  degree  of  Pi{x),  or  that  of 
some  of  its  neighbours,  whenever  an  edge-detection  mechanism  finds  an  edge  at  the  ith 
interval.  By  lowering  the  degree  of  some  polynomial  pieces  close  to  an  edge,  one  can 
avoid  crossing  the  edge  in  the  prediction  step,  as  much  as  possible.  This  option  is  closely 
related  to  the  nonlinear  multiscale  transformation  considered  in  [6]  (within  the  Lifting 
framework),  where  the  nonlinearity  comes  in  from  adaptively  choosing  from  the  (1  ,M) 
family  of  linear  filters. 


Once  Si  is  determined  (i  € Si),  Pi(x)  can  be  uniquely  determined  when  degree  Pi(x)  = 
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#Si  [1]  so  that 

^ l*xin  1 

7 / pi(x)dx  — v^1  for  me  Si,  (2.11) 

hJ-l  JaJ-\ 

and  the  prediction  operator  is  then  defined  by  (2.10). 

One  can  be  slightly  more  ‘sophisticated’  in  the  design  of  the  polynomial  pieces.  The 
Subcell  Resolution  technique  [4,  13]  allows  to  account  for  discontinuities  within  a cell  as 
follows.  If  an  edge  is  detected  in  the  ith  cell,  the  polynomial  piece  p,  (x)  is  discarded  and 
substituted  by  its  left  and  right  neighbours,  pi+ \(x)  and  p;_ i(rr),  assuming  that  their 
respective  stencils  do  not  intersect,  i.e.  iS;_i  H^+i  = 0-  At  a true  one-dimensional  edge 
(a  jump)  on  the  ith  cell,  the  function 

1 fv  1 rx2< 

F(y)  = -r-  pi-i{x)dx+—\  pi+1(x)dx 

”■]  Jx32i_  J "j  Jy 


will  have  a zero  on  the  ith  cell  [13],  say  rj,  and  the  location  of  rj  is  used  to  substitute  the 
polynomial  piece  Pi{x)  by  the  discontinuous  piecewise  polynomial  function 


Pi-i  0*0 
Pi+i(x) 


X<T), 
X > Tj. 


(2.12) 


The  prediction  operator  is  again  defined  by  (2.10)  at  nonsingular  cells  (cells  in  which  no 
edge  has  been  detected),  while  at  the  singular  cell 


■ i 1 fx™  • , 1 fx 

(PjV3  ) 2i  = 7-  qi(x)dx,  ( PjV 3 I)2i-1  = r-  5»(.' 

h3  Jxi2i_1  hi  Jx’2i_2 


(x)dx. 


In  practice  it  is  unnecessary  to  compute  explicitly  the  value  of  r)\  only  its  location  with 
respect  with  *2  i-  .j  is  needed,  which  can  be  found  by  a sign  check.  We  refer  the  reader 
to  [4]  (and  references  therein)  for  specific  details  on  this  technique,  in  particular  on  the 
detection  mechanism,  and  on  its  performance. 

3 The  question  of  stability:  Error  control  versus  synchroniza- 
tion, with  numerical  examples 

Lossy  coding  schemes  introduce  errors  into  the  transform  coefficients,  and  it  becomes 
crucial  that  the  nonlinearities  do  not  unduly  amplify  these  errors.  In  lossy  compression 
the  decoder  only  has  the  quantized  detail  coefficients.  If  we  use  a nonlinear  prediction 
operator  (whether  it  is  constructed  as  described  in  the  previous  section  or  based  on 
locally  adapted  filters,  as  in  [6]  within  the  Lifting  framework),  the  quantization  errors  in 
coarse  scales  could  cascade  across  the  scale  ladder  and  cause  a series  of  incorrect  choices 
(either  on  the  filters  or  on  the  stencils)  leading  to  serious  reconstruction  errors. 

To  avoid  incorrect  choices  in  the  prediction  step,  whether  within  Harten’s  or  the 
Lifting  framework,  one  would  need  to  send  side  information  on  which  filter  was  used 
(Lifting)  or  what  was  the  interpolatory  stencil  (Harten’s).  This  is  clearly  inappropriate 
when  trying  to  design  a compression  scheme.  One  way  to  avoid  storing  (and  sending)  side 
information  is  to  somehow  synchronize  the  nonlinear  prediction  operators  in  the  encoder 
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and  the  decoder,  so  as  to  ensure  that  at  a given  spatial  location  on  a given  scale,  the 
prediction  operator  will  select  the  same  stencil  (filter  bank),  both  in  the  encoding  and 
the  decoding  steps. 

Within  the  Lifting  framework,  synchronization  is  achieved  in  [6]  by  changing  the 
typical  Split-Predict-Update  steps  to  Split-Update-Predict.  In  doing  so,  it  is  possible  to 
base  the  choice  of  predictor  directly  on  already  ‘quantized  data’,  thus  synchronizing  the 
nonlinear  decisions  made  by  the  encoder  and  the  decoder. 

Within  Harten’s  framework,  synchronization  is  just  a consequence  of  a strategy  that  is 
designed  to  fully  control  the  compression  error.  Because  the  main  design  tool  in  Harten’s 
framework  for  multiresolution  is  a reconstruction  technique,  and  because  A.  Harten  had 
already  worked  with  nonlinear  reconstruction  techniques  in  the  context  of  the  numerical 
simulation  for  hyperbolic  conservation  laws,  so-called  Error- Control  (EC)  strategies  can 
be  found  already  in  the  early  papers  of  Harten  on  multiresolution  [14]. 

Harten’s  mechanism  to  control  the  global  accumulated  error  is  based  on  a modification 
of  the  direct  multiscale  transformation,  M,  that  ensures  a prescribed  tolerance  on  the 
global  prediction  errors  (explicit  error  bounds  can  be  found  in  [4,  13]).  The  modified 
transformation  incorporates  the  quantizer  to  the  direct  multiscale  transformation  in 
such  a way  that  the  prediction  operator  in  the  encoder  also  acts  on  already  ‘quantized’ 
data,  hence  synchronization  is  achieved  because  the  nonlinear  prediction  operators  both 
in  M and  M_1  work  on  the  same  set  of  discrete  data  at  each  resolution  level. 

To  illustrate  the  effect  of  the  different  techniques,  we  take  a particular  nonlinear 
prediction  operator,  a third  order  ENO  reconstruction  technique  with  Subcell  Resol- 
ution, as  described  in  last  section.  We  denote  by  Msn  the  multiscale  transformation 
(2.6),  while  Mgf t denotes  the  EC  modified  transform  as  described  in  [2,  4],  and  M§R  a 
multiscale  transformation  in  which  only  synchronization  is  enforced,  as  proposed  in  [6] . 
The  quantization  step  is  carried  out  as  follows: 

qu(tf')  = 2e,-round  \d? /(2ej)\ 

and  it  is  incorporated  to  the  direct  transformation  in  M™R  and  M §R  (see  [2,  6]  for 
specific  details),  while  in  MRr  it  is  applied  to  the  scale  coefficients  obtained  after  the 
transformation.  In  the  numerical  tests  we  report,  we  take  cl  = 8 with  L = 4 and 
£j  = £j+i /2- 

We  consider  two  different  images:  the  familiar  image  of  Lena  as  an  example  of  a ‘real’ 
image,  and  a purely  geometrical  image,  to  which  texture  has  been  added,  as  in  [6]. 

After  the  direct  transformation  (plus  the  quantization  step)  has  taken  place,  a lossless 
Lempel-Ziv  compression  algorithm  is  applied  to  reduce  the  size  of  the  transformed  image, 
then  a compression  ratio  is  computed  as  the  number  of  bits  of  the  compressed  repres- 
entation over  the  number  of  bits  of  the  original  image.  To  recover  the  original  image,  we 
undo  the  lossless  compression  and  transform  back  using  (2.7)  in  all  three  cases.  The  full 
compression  algorithm  is  identified  in  each  case  by  an  acronym,  ‘ST’  for  Msr , ‘EC’  for 
and  ‘SYNC’  for  M§R. 

In  Tables  1 and  2 we  compile  a number  of  quantities  that  measure  the  ‘quality’  of  the 
reconstructed  image,  and  therefore  the  robustness  and  reliability  of  each  multiresolution- 
based  compression  algorithm,  the  magnitude  of  the  global  compression  error,  measured 
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Method 

II  -Hoc 

II -lb 

II  ’ 1 1 2 

rc 

entropy 

ST 

258 

5.71 

9.08 

11.3:1 

.6449 

SYNC 

195 

6.45 

9.82 

7.9:1 

.8875 

EC 

25.4 

4.47 

5.73 

9.7:1 

.6850 

Tab.  1.  Geometrical  image. 


Fig.  1.  Geometrical  image:  (a)  original,  (b)  ST,  (c)  EC,  (d)  SYNC. 


in  various  norms,  the  compression  rate  rc  and  the  entropy  of  the  transformed  image. 
The  reconstructed  images  in  both  cases  can  be  observed  in  Figures  1 and  2. 

It  can  be  clearly  observed  that  the  absence  of  any  type  of  synchronization  procedure 
can  lead  to  a very  poor  reconstructed  image.  Synchronization  only  improves  the  quality, 
but  is  not  as  robust  as  the  full  EC  mechanism,  designed  in  this  case  to  enforce  a certain 
error  bound  in  the  2-norm  (as  observed  in  Tables  1 and  2,  the  2-norm  of  the  global  error 
is  kept  below  ez,  = 8).  It  is  worth  mentioning  that  the  compression  rate  and  the  entropy 
of  the  compressed  data  are  all  very  dose,  however  the  visual  quality  of  the  reconstructed 
image  is  significantly  better  for  the  EC  compression  algorithm. 
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